#ifndef _LSTATFIT_H_
#define _LSTATFIT_H_

#include <clx/libECM.h>

#include <gsl/gsl_fit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

struct _msl_fit_t;

typedef ECMBOOL (*data_reader_t)(char *buf, ecm_double *pxvalue, ecm_double *pyvalue);

typedef int (*dist_func_f_t)(const gsl_vector *x, void *data, gsl_vector *f);
typedef int (*dist_func_df_t)(const gsl_vector *x, void *data, gsl_matrix *J);

typedef void (*dist_report_func_t)(FILE *fp, struct _msl_fit_t *fit, ecm_uint32 xdiff);
typedef ecm_double (*dist_goodness_func_t)(struct _msl_fit_t *fit);

typedef enum {
	MSL_FIT_NONE = 0,
	MSL_FIT_EXPONENTIAL = 1,
	MSL_FIT_BETA,
	MSL_FIT_LOGNORMAL,
	MSL_FIT_WEIBULL,
	MSL_FIT_MAX = MSL_FIT_WEIBULL
} msl_fit_type_t;

typedef struct {
	const char	*desc;
	dist_func_f_t	dist_f;
	dist_func_df_t	dist_df;
	dist_report_func_t	dist_report;
	dist_goodness_func_t	dist_goodness;
	ecm_uint32	p;
	ecm_double	x_init[5];
	ecm_double	tolerance;
} fit_info_t;

typedef struct _msl_fit_t {
	msl_fit_type_t	fit_type;
	fit_info_t	*info;
	char	*name;
	const gsl_multifit_fdfsolver_type	*T;
	gsl_multifit_fdfsolver	*s;
	gsl_multifit_function_fdf	f;
	gsl_matrix	*covar;
	ecm_uint32	n;
	ecm_double	xrange;
	ecm_double	max;
	ecm_double	*xvalues;
	ecm_double	*yvalues;
	unsigned int	iter;
	unsigned int	status;
} msl_fit_t;

extern ecm_uint32	n_inputs;
extern ecm_double	*xvalues;
extern ecm_double	*yvalues;

ECMBOOL read_input(const char *inputfile, data_reader_t reader);

extern msl_fit_t *lstat_fitting(const char *name,
				msl_fit_type_t fit_type,
				ecm_uint32 n_inputs,
				ecm_double *xvalues,
				ecm_double *yvalues);
void lstat_fitting_free(msl_fit_t *fit);

void lstat_report_fitting(msl_fit_t *fit, ecm_uint32 xdiff);
ecm_double lstat_get_goodness(msl_fit_t *fit);

msl_fit_type_t parse_model_type(const char *typestr);

#endif
